import { IPoint } from "@/models/gui";
import { IFn } from "@/models/math";
import { getRidgeRegressionFittingFn } from "./ridge-regression";

// 多项式插值 手动求解
export function getInterpolationFn (points: IPoint[]): IFn {
  const n = points.length;
  const B: number[][] = [];
  const C: number[] = points.map(() => null);
  points.forEach((p, i) => {
    B[i] = [1];
    for (let j = 1; j < n; j++) {
      B[i][j] = B[i][j - 1] * p.x;
    }
  });
  const Y = points.map(p => p.y);
  for (let k = 0; k < n - 1; k++) {
    const bk = B[k];
    for (let i = k + 1; i < n; i++) {
      const bi = B[i];
      const b = bi[k];
      for (let j = k; j < n; j++) {
        bi[j] = bi[j] / b - bk[j] / bk[k];
      }
      Y[i] = Y[i] / b - Y[k] / bk[k];
    }
  }
  C[n - 1] = Y[n-1] / B[n-1][n-1];
  for (let k = n - 1; k >= 0; k--) {
    const bk = B[k];
    let sum = 0;
    for (let j = k + 1; j < n; j++) {
      sum += bk[j] * C[j];
    }
    C[k] = (Y[k] - sum) / bk[k];
  }
  return (x) => {
    let y = 0;
    let xn = 1;
    C.forEach(c => {
      y += c * xn;
      xn = xn * x;
    });
    return y;
  };
}

// 插值 (x - x0)(x - x1)...形式
export function getInterpolationFn2 (points: IPoint[]): IFn {
  const X = points.map(p => p.x);
  const Y = points.map(p => p.y);
  const B = [];
  const n = points.length;
  for (let i = 0; i < n; i++) {
    B[i] = Y[i];
    const xi = X[i];
    for (let j = 0; j < n; j++) {
      if (j != i) {
        B[i] /= xi - X[j];
      }
    }
  }
  return (x) => {
    const index = X.indexOf(x);
    if (index >= 0) {
      return Y[index];
    }
    let A = 0;
    let C = 1;
    X.forEach((xi, i) => {
      C *= x - xi;
      A += B[i] / (x - xi);
    });
    return A * C;
  };
}

// 回归分析 手动求解
export function getCurveFittingFn1 (points: IPoint[]): IFn {
  function h (f1: IFn, f2: IFn, X: number[]): number {
    let result = 0;
    X.forEach(x => {
      result += f1(x) * f2(x);
    });
    return result;
  }
  const n = points.length - 1;
  const X = points.map(p => p.x);
  const Y = points.map(p => p.y);
  const A: number[] = [];
  const cores: IFn[] = [];
  for (let i = 0; i < n; i++) {
    cores[i] = (x) => Math.pow(x, i);
  }
  const f: IFn = (x) => {
    const index = X.indexOf(x);
    return Y[index] || 0;
  };
  const B: number[][] = [];
  for (let i = 0; i < n; i++) {
    B[i] = [];
    for (let j = 0; j < n; j++) {
      B[i][j] = h(cores[i], cores[j], X);
    }
    // console.log([...B[i]]);
  }
  const C: number[] = [];
  for (let i = 0; i < n; i++) {
    C[i] = h(f, cores[i], X);
  }

  // console.log([ ...C ]);
  for (let k = 0; k < n - 1; k++) {
    const bk = B[k];
    for (let i = k + 1; i < n; i++) {
      const bi = B[i];
      const b = bi[k];
      for (let j = k; j < n; j++) {
        bi[j] = bi[j] / b - bk[j] / bk[k];
      }
      C[i] = C[i] / b - C[k] / bk[k];
      // console.log([ ...bi ]);
    }
  }
  // console.log([ ... C]);
  A[n - 1] = A[n-1] / B[n-1][n-1];
  for (let k = n - 1; k >= 0; k--) {
    const bk = B[k];
    let sum = 0;
    for (let j = k + 1; j < n; j++) {
      sum += bk[j] * A[j];
    }
    A[k] = (C[k] - sum) / bk[k];
  }
  // console.log(C);
  return (x) => {
    let y = 0;
    A.forEach((a, i) => {
      y += a * cores[i](x);
    });
    // console.log(x, y);
    return y;
  };
}

// 矩阵运算库 岭回归
export function getCurveFittingFn2 (points: IPoint[]): IFn {
  const cores: IFn[] = [];
  for (let i = 0; i < points.length - 1; i++) {
    cores[i] = (x) => Math.pow(x, i);
  }
  return getRidgeRegressionFittingFn(cores, points);
}